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Abstract 

Structure formation within the Lemaitre-Tolman model is investigated in a gen- 
eral manner. We seek models such that the initial density perturbation within a 
homogeneous background has a smaller mass than the structure into which it will 
develop, and the perturbation then accretes more mass during evolution. This is a 
generalisation of the approach taken by Bonnor in 1956. It is proved that any two 
spherically symmetric density profiles specified on any two constant time slices can 
be joined by a Lemaitre-Tolman evolution, and exact implicit formulae for the arbi- 
trary functions that determine the resulting L-T model are obtained. Examples of 
the process are investigated numerically. 
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1 Introduction 

Though the Lemaitre-Tolman (L-T) model has been studied extensively (see the ques- 
tion of whether galaxies and other modern cosmic structures can grow from small initial 
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perturbations, using this exact inhomogenous cosmological model, has not been clearly 
answered. 

Long ago Bonnor P| considered a version of the general problem, in which the pertur- 
bation consisted of an interior region matched to a Friedmann-Lemaitre- Robertson- Walker 
(Friedmann) exterior, so the initial fluctuation of density included all the dust particles 
that would enter the future galaxy, and so the outer edge of the perturbation had to be 
comoving with the (spatially homogeneous) background flow ever after. Because of this, 
and given the present age of the Universe, the initial fluctuation had to have an ampli- 
tude many times larger than a statistical fluctuation could have. However current matter 
models, that allow perturbations to grow before recombination, have successfully predicted 
temperature perturbations in the CMB of order 10 -5 . 

The very existence of inhomogeneous cosmological models (i.e. spatially inhomogeneous 
solutions of Einstein's equations with expanding matter), such as the L-T |3|, |J or Szekeres 
II U models, shows that non-Friedmannian distributions of density and velocity would 
have been coded in the Big Bang and need not be "explained" as statistical fluctuations 
that appeared within a homogeneous background during evolution. Moreover, since the L- 
T collection of models is labelled by two arbitrary functions of mass, that reduce to specific 
forms in the Friedmann limit, it follows that the dust Friedmann models are a subset of 
measure zero within the L-T set. Consequently, the Friedmann models are very improbable 
statistically and, assuming that our physical Universe is homogeneous indeed, one needs 
to explain how homogeneity might have come about out of inhomogeneous initial data, 
not the other way round. This is what inflation is believed to have done. However, in this 
paper we shall accept a high degree of homogeneity at decoupling, and we will determine 
how fast condensations can grow, once they appear in a homogeneous background. 

In this paper we shall relax Bonnor's assumptions, and consider the general case. In 
particular, we envisage a scenario in which the mass of the initial fluctuation is much 
smaller than the mass of the condensation into which it will develop, and that it captures 
more mass during its evolution.^] The outer edge of the growing condensation will thus 
not be comoving with the background flow. The calculations, although based on exact 
formulae, will have to be carried out numerically. 

In the following we set A = 0, as its effect is felt primarily at late times over the long 
range, and so will not strongly affect structure formation. Also, current interpretations 
of the CMB and Supernova data that estimate a non-zero Q^, should be regarded as 
provisional, since several reasonable alternatives have been put forward. 



1 Such a modification of Bonnor's method was suggested by S. Bazariski during one of the seminars by 
A. K. in Warsaw. 
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2 Basic properties of the Lemaitre-Tolman model. 



The Lemaitre-Tolman (L-T) model || is a spherically symmetric nonstatic solution of 
the Einstein equations with a dust source. Its metric is: 

ds 2 = dt 2 ^zrdr 2 - R 2 (t,r)(d$ 2 + sin 2 ^ 2 ), (2.1) 

1 + 2E{r) 

where E(r) is an arbitrary function (arising as an integration constant from the Einstein 
equations), R, r is the derivative of the function R(t,r) by r, and R obeys the equation 

R,t 2 = 2E + 2M/R+ -AR 2 , (2.2) 

3 



where A is the cosmological constant. Eq. Q2.2f) is a first integral of one of the Einstein 
equations, and M(r) is another arbitrary function that arises as integration constant. The 
matter-density is: 

K ,= |^, where K = (2.3) 



R 2 R' c 



In the following, we will assume A = 0. Then eq. (|2.2|) can be solved explicitly. The 
solutions are: 



When E < 0: 



. . M . . 
R\t,r) = -— (l-cosr?), 

V - sin 77 = ( ~ 2 ^ )3/ (t - t B (r)). (2.4) 



where rj is a parameter; when = 0: 



i?(i,r) 



^M(t-t B (r)) 2 



1/3 



(2.5) 



and when E > 0: 



M 

= ^(coshr] - 1), 

sinhr; - 77 = V } M (t - t B (r)), (2.6) 

where te(r) is one more arbitrary function (the bang time). Note that all the formulae 
given so far are covariant under arbitrary coordinate transformations r = g(r'), and so r 
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can be chosen at will. This means one of the three functions E(r), M(r) and £b(V) can be 
fixed at our convenience by the appropriate choice of g. 

In a general L-T model, E may change sign, having both re-collapsing and ever- 
expanding regions. Also the space t = const with E(r) < everywhere is not necessarily 
closed (and the one with E(r) > is not necessarily infinite), see Refs. [0, §. 

The Friedmann models are contained in the Lemaitre-Tolman class as the limit: 



t B = const, \E\ 3/2 /M = const, 



(2.7) 



and one of the standard radial coordinates for the Friedmann model results if the coordi- 
nates in (p.4|) - (|2.6|) are chosen so that: 



M = M r 3 , (2.8) 

where Mq is an arbitrary constant; so that E = Eqt 2 , Eq = const. 

It will be convenient in most of what follows to use M(r) as the radial coordinate (i.e. 
r' = M(r)) because in the structure formation context one does not expect any "necks" or 
"bellies" where M, r = 0, and so M(r) should be a strictly growing function in the whole 
region under consideration. (See the papers by Barnes and, especially, by Hellaby || [10 



for descriptions of necks. Some properties of a neck appeared also in the paper by Novikov 
TJ). Then: 



k P = 2/(R 2 R, m ) =6/(i? 3 ), M . (2.9) 

We also, in searching for realistic models, prefer L-T models that are free of shell 
crossings M. This is not because shells of matter cannot collide, but because the L-T 
co-moving description breaks down there. 

We shall now apply this model to the problem of structure formation within the exact 
relativity theory (i.e. without approximations). We believe this question has not been 
satisfactorily answered so far, and so it deserves to be investigated more thoroughly, both 
as an important consequence of L-T models, and with a view to possible cosmological 
applications. We seek accreting models in which a small initial fluctuation at decoupling 
captures more mass during its evolution, thus growing in extent as well as in density 
contrast. 
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3 The evolution as a mapping from an initial density 
to a final density. 



The evolution of the L-T model is usually specified by defining the initial conditions — the 
distributions of the Big Bang time tg(M) and of energy E(M), or by specifying e.g. the 
density p(£i, R) and velocity R t t(ti, R) at an initial instant t — t±. It is, however, possible, 
to approach the problem in a different way: to specify the density distributions at two 
different instants, t — t\ and t = t 2 , calculate the corresponding E(M) and te(M), and in 
this way obtain a definite model. It is not immediately obvious whether all pairs of density 
distributions may be connected by an L-T evolution of a chosen type; nor whether one can 
ensure shell-crossings do not occur between t\ & t 2 . However, such a mapping from an 
initial density to a given final density should exist in many cases, especially with a sensible 
choice of density profiles. 

In fact, it will be proven below that any initial value of density at a specific position 
(r, M const) can be connected to any final value of density at the same position by one 
of the Lemaitre-Tolman evolutions (either E > 0, or E < 0, or, in an exceptional case, 
E — 0). In the Friedmann limit, any two constant densities can be connected by one of 
the fc>0,A;<0or&; = Friedmann evolutions. 

For definiteness, it will be assumed in the following that the final instant t 2 is later than 
the initial instant t\, i.e. t% > t%, and that the final density p(t2,M) is smaller than the 
initial density p(t±,M) at the same M. We thus assume that matter has expanded along 
every world-line, but the proof can be easily adapted to the collapse situation. 

3.1 Hyperbolic regions 

Let us consider the L-T model with E > 0. Let the initial and final density distributions 
at t = t\ and t = t 2 be given by: 

p(t u M) = Pl (M), p{t 2 , M) = p 2 {M) (3.1) 

From (|2.3|) we then have, for each of t± Sz t^- 

r M ft 

R 3 (U, M) - R 3 mini = / -— — dM' := R t 3 (M), . = 1,2 (3.2) 

JM min Kpi(M') 

and Ri(M) > R\(M) in consequence of p(t 2 ,M) < p(ti,M). In the following we will 
assume there is an origin where M = and R(ti,0) = 0, so that i? mini = = M min is 
validf]. Solving the two parts of Q2.6|) for t(R, r) and writing it out for each of (ti, R\) and 

2 For examples where this is not the case, see the papers by Einstein & Straus JlJ, with M min ^ 
& R m i n — 0, and by Bonnor & Chamorro H] where M = from R = to R m . m . Also "necks" are the 
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(fa, R2) leads to: 



M 

ts — fa 



1 + 2ERi/M) 2 - 1 - arcosh(l + 2ERi/M) 



(2Efl 2 

and then eliminating t B between the two versions of (Of) we find: 



i = 1,2 (3.3) 



1 + 2ER 2 /M) 2 - 1 - arcosh(l + 2ER 2 /M) 
1 + 2ER 1 /M) 2 - 1 + arcosh(l + 2ER 1 /M) = [(2E) 3/2 /M](fa - fa), (3.4) 

We shall prove that this equation has one and only one solution E(M) > (in addition to 
the trivial solution E = 0) provided that fa and fa obey a certain inequality (see below). 
In fact, the inequality will exclude the E < models. 

For ease of calculations, let us denote: 

x := 2E/M 2/3 , cii = Ri/M 1/3 , z = 1, 2; 



^ff(x) := y (1 + a 2 x) 2 — 1 — arcosh(l + a 2 x) — y (1 + a x x) 2 — 1 + arcosh(l + dx) 

-(fa ~ fa)x 3/2 := xn(a:) - (fa - fa)x 3/2 , (3.5) 

Our problem is then equivalent to the following question: for what values of the parameters 
(12 > ai and fa > fa, does the equation ipn(x) = have a solution 1 ^ 0? Note that ipn 
has a zero at x = that the more correct but less convenient x~ 3 ^ 2 Xh — (fa — fa) does not 
have. 

Since x > and > by definition, iPh(x) and XffO^) are well-defined for any 
.x G [0, 00). Note that 

, arcosh(l + a^x) arcosh(l + a>ix) , . 

hm . = lim — q/Q = 0, (3.6) 

x ^°° J (I + a, t x) 2 - 1 x - >0 ° (t 2 -ti)x 3 / 2 



and so in determining the sign of iPh(x) as x — * 00, the arcosh-terms can be neglected. 
Note also 

lim . XH ^\ /9 = 0. (3.7) 
s-00 (t 2 - fa)x 3 / 2 v ; 

Hence, the last term in ipH becomes dominant when x — > +00, and so 

lim iPh(x) = —00. (3.8) 



locus of a minimum in M and R. 
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It is easy to see that i[>h(0) = 0, but we wish to know whether ipH — > + or ipu ^ as 
x — > 0. For this purpose, note that 



-^-ip H (x) = ^H,x 
ax 



«2 



3/2 



3/2 



; ; - -(t 2 ~ h) 



(3.9) 



from which it follows that 



^h,x(0) = 0, lim = -oo. (3.10) 



It is also easy to see that the term in square brackets in ( |3.9| ) is a strictly decreasing 
function for all x G [0, oo), and so it may be equal to zero in at most one point. Since 
it goes to the negative value — (3/2) (£2 — ti) when x —>■ 00, it will have a zero when it is 
positive at x = 0, i.e. when 

h - h < ^ (a 2 3 / 2 - a^ 2 ) . (3.11) 

By comparison with (|2.5|) , since Oj = Ri/M 1 ^ 3 , the above is seen to be equivalent to the 
statement that between t\ and t 2 , R(t,M) has increased by more than it would have 
increased in the E = L-T model. This is a necessary condition for the existence of an 
E > evolution connecting R(t X) M) to R(t 2 ,M). 



It is also a sufficient condition, as we now explain. With ( |3.11| ) fulfilled, lim^o ^h, 



+ , i.e. ipH,x > in a neighbourhood of x = 0, then it goes through zero exactly once, at 
some x = x m , and becomes negative. This means that iPh{x) itself is increasing from the 
value at x = 0, to a maximum at x = x m , and is then decreasing all the way to x — > 00 
where it becomes —00. Hence, at one and only one x = x > x m , i/jh(xo) — 0. This implies 
that eq. ( p.4|) defines a function E(M) in the whole range of M in which (|3.11|) is fulfilled. 
Examples of the functions iPh(x) from (3.8) that obey or do not obey ( 3.11|) are shown in 
Fig. 1. 



Fig. 1. goes here 



There remains a practical problem for the numerical calculation of E(M). Since the 
range of x is infinite in (|3.5|) , an initial value xa < 00 such that i/)h{%a) < has to be 
determined first. For this purpose, note that for large x we have \J (I + dix) 2 — 1 ~ l + djX. 
Together with (|3T6| ), this implies that iPh(x) is well approximated for large x by 

i/; A (x) := (a 2 - a x )x - {t 2 - t^x^ 2 . (3.12) 
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Indeed, it is easy to verify that for all x > 



ip H {x)<ijj A {x). (3.13) 

[Writing ip H (x) - ip A (x) = Xh{x) - (a 2 - a x )x, and (ip HjX - ip A>x ) = W(a 2 ) - W(a x ), 
where W(a) = a 3 ^ 2 y/x/^2 + ax — a, we note that W(a) is a decreasing function of a for 
every x > 0, while ipni^) = V'a(O) = 0. Hence, iPh{.x) < i/j a (x) for x > 0.] Therefore, if 
?AaOa) = 0, then ^ H {x A ) < 0. 

The solution of ip A (x) = is 

and so a;^ is a good initial value for the numerical program that will find a solution of 
^h{.x) = by bisecting the segment [0,xa] and checking the sign of i/jh(x a /2). 

Also, for numerical purposes, the limits of some of the functions at M = must be 
calculated separately, as explained in appendix [TJ. 

3.2 Still-expanding elliptic regions 

For E < 0, a similar result holds, but with one more refinement: depending on the value of 
(^2 — the final density will be either in the expansion phase or in the recollapse phase 
(and only in one of these phases). The dividing value of (t 2 — ti) will come out in the proof 
below. 

Let us assume that the i] of (|2.4|) is in [0,7r] for both values of £«, so that the final 
density is still in the expansion phase of its evolution. (For rj G [k, 27r], the solutions for 
(ti — t B ) are different, and they will be considered separately below.) The analogs of eqs. 
<K1) and (ET3) are then: 



M 

tB — ti ~ 



[-2Ef/ 2 
and 



arccos(l + 2ERi/M) -Jl-(1 + 2ER { /M) 



(3.15) 



M*) = °, ( 3 - 16 ) 

where this time 



x 



-2E/M 2/3 , 



ipx( x ) '■= arccos(l — a 2 x) — Jl — (1 — a 2 x) 2 — arccos(l — a±x) + Jl — (1 — a\x) 2 
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-(*2 - h)x 3/2 := - (t 2 - ti)x 3/2 , 



(3.17) 



the definitions of being still ( |3.5| ). 

The reasoning is entirely analogous to the one for (|3.5|) , but this time the arguments 
of arccos must have absolute values not greater than 1. This implies x < for both i, 
and so, since a 2 > a\ 



< x < 2/a 2 , 



(3.18) 



which means: if there is any solution of ( |3.16| ), then it will have the property ( |3 . 1 8[ ) . The 
two square roots in ( |3.17| ) will then also exist. Eq. ( |3.18| ) is equivalent to the requirement 
that (R, t ) 2 (in Q2.2D with A = 0) is nonnegative at both t\ and t 2 . 

Note that 



^x(O) = 0, 



■0x(2/o2) —7t — arccos(l - 2a 1 /a 2 ) + 2^Ja 1 /a 2 - (ai/a 2 ) 2 - (2/a 2 ) 3/2 (t 2 - ti), 



ax 



«2 



3/2 



a i 



3/2 



V2 - a 2 x \/2~ 



3, 



(3.19) 
(3.20) 



In consequence of a 2 > a%, the term in square brackets is now an increasing function of 
x, and it becomes +oo at x = 2/a 2 (which only means that ipx{x) has a vertical tangent 
there). Hence, the term can go through zero at most once, and it will do so when it is 
negative at x — 0, i.e. when the opposite to (|3.11|) holds 



V2 ( 
3 



h - U > — (a 2 3/2 



a\ 



3/2 



(3-21) 



This means that the model must have expanded between ti and t 2 by less than the E = 
model would have done. If (|3.21|) does not hold, then ipx,x is positive for all x > 0, which 
means that ipx{x) is increasing and will not be zero for any x > 0. Hence, ( |3.21|) is a 



necessary condition for the existence of a solution of (|3.16|) . 

With (|3.21| ) fulfilled, ipx,x{x) becomes negative for some x > 0, then goes through zero 
exactly once and then is positive all the way up to x = 2/a 2 . This implies that ipx(x) 
initially decreases below 0, then has exactly one minimum and is increasing up to the value 
( |3.19| ) at x = 2/a 2 . Hence, ipx{x) will have a zero for x > only if ipx{2/a>2) > 0, i.e. if 



t 2 -h< {a 2 /2 f 2 



71 



arccos(l — 2a\ja 2 ) + 2Ja\ja 2 — (ai/a 2 ) 



(3.22) 



The inequality ( |3.22|) is consistent with (|3.21|) , see appendix ^. Eqs. (|3.21|) and ( |3.22|) 
together are a necessary and sufficient condition for ( |3.16| ) - (|3.17|) to define a function 
E(M) < for which R(t 2 , M) is still in the expansion phase of the model. 
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3.3 Recollapsing elliptic regions 



The reasoning above applied only in the increasing branch of R in 
branch, where 77 G [n, 2n], instead of Q3.16j ) - (|3.17|) we obtain 



For the decreasing 



ts — t\ — 
= h- 



M 



_ 2 £)3/2 

M 
'-2Efl 2 



arccos(l + 2ER X JM) - Jl - (1 + 2ER 1 /M) 



tt + arccos(-l - 2ER 2 /M) + Jl - (1 + 2ER 2 /M) 2 



(3.23) 



and 



tpc — 0, where 



ipc(x) '■= tt + arccos(— 1 + a 2 x) + yl — (1 — a 2 x) 2 — arccos(l — a\x) + yl — (1 — a x x) 2 

-(t 2 - tl )xV 2 . (3.24) 



The derivative of this is 



d_ 

dx 



lp c (x) = tpc, x = -yfx 



«2 



.3/2 



+ 



3/2 



y/2 — a 2 x a/2 — a\X 2 



(3.25) 



and is negative for all x > 0. Since 



Vc(0) = 2tt > 0, 



(3.26) 



the solution of ipc{ x ) — for x > will exist if and only if ipci^/ a 2) — 0, which translates 
into the opposite of (|3.22|) : 



*a - *i > (« 2 /2) 



3/2 



7r — arccos(l — 2a\ja 2 ) + 2\ja\ja 2 — (ai/a 2 ) 2 



(3.27) 



Thus the two densities can be connected by an E < L-T evolution that is recollapsing at 
time t 2 if ( ET27D is obeyed. 

Examples of functions ipx(x) obeying or not obeying Q3.21Q and (|3.22 ), and of functions 
ipc( x ) obeying or not obeying ([3.271) are shown in Figs. 2 and 3. 



Figs. 2 and 3 go here 
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3.4 Summary 



The above analysis considered only single world-lines, that is, single M values. We extend 
this to the whole of Pi{M) by noting that E(M) and £b(M) are arbitrary functions in the 
L-T model, and so continuous pi will generate continuous E & £g. 

The meaning of the limiting cases is now easy to understand. In ( 3.11|) , at M values 
where t 2 — £1 = (v^2/3)(a2 3//2 — ai 3 ^ 2 ), the final state results from the initial one by a 
parabolic (E = 0) evolution, so this M value is on the boundary between an elliptic region 
and a hyperblic one. Eq. Q follows as the E -> limit of O and of In ( gjgfl 

and ( |3.27| ), for M values where the equality holds, the final state is exactly at the local 
moment of maximal expansion, separating a region of P2 that is already recollapsing from 
one that is still expanding. 

However, when E < it must be remembered that the signature of the metric requires 
that 



E(M) > -1/2, (3.28) 



and so, once E(M) has been calculated, (|3.28|) will have to be checked. Note that the 



Friedmann model in standard coordinates has exactly this problem — with IE = —kr 2 
and M = M r 3 , blindly continuing through r = 1 will make E < —1/2 and M > Mjj niverse . 
Indeed, given two uniform densities P2 < pi that are appropriate for a closed Friedmann 
model, the integral fl3.2p can be extended to arbitrarily large M & R. Thus the occurrence 
of E = —1/2 is not a problem, but rather an indication that the maximum of the spatial 
section has been reached. One should record the values of R m ax and M max , and then use 

3L» - R 3 (U, M) = / -— — dM'. (3.29) 

and distinguish the M values beyond the maximum from those in front of it. 

Another serious possibility is that shell crossings, where the density diverges and 
changes sign, may occur. If they occur between t\ & t 2 the model evolution is unsat- 
isfactory, but if they occur before £1 or after t 2; this may not be of much concern. The 
conditions on E(M) & Ib{M) for avoiding them [[7] must also be checked. 

All these considerations apply to the Friedmann limit, but it must be remembered that 
in comparing models in a continuous Friedmann family, one must not scale the curvature 
index A; to +1 or —1 when it is nonzero. The parameter k is adapted to the initial and 
final densities together with Mq. With k scaled to ±1, taking the limit k — > within the 
family becomes impossible, and the inequalities ( ft.ll ) and (|3.21|) do not come up. 



In summary, for densities p 2 {M) < pi(M) at times t 2 > £1 we have a 2 > a± where 
Oj = i?j/M 1//3 , and writing 

a = d\ja 2 
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the nature of the LT model that evolves between these states at a given M is: 



Hyperbolic E > : 
If 

t 2 -h < (V2a 3 2 /2 /3)(l-a^ 2 ) 

then 

E = xM 2/3 /2 

where x solves 

= iPh(x) = \J {1 + a 2 x) 2 — 1 — arcosh(l + a 2 x) 

-J {I + aix) 2 - 1 + arcosh(l + a x x) - (t 2 - t^x^ 2 



and 



3/2 



y^l + ajx) 2 — 1 — arcosh(l + a,ix) 



Parabolic E = : 
If (t 2 — ti) is close to 

t 2 -t 1 = {y/2a\ 12 /3)(l-a 3 / 2 ) 
then a series expansion gives 

£ = x M 2/3 /2 

where x solves 

. . V2 X 3 / 2 f 3/2 / 3 9 2 2 N 

y w 3 I 2 V 20 224 2 , 



and 



v^.3/2/. 3 9 2 2 



t B m h a/ 1 <2j£ H fyx' 

3 1 V 20 224 ' 
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Elliptic E < and still expanding at U: 
If 




(a 2 /2) 3/2 [7r - arccos(l - 2a) + 2y/a - a 2 ] > t 2 - h > (V2 a 2 /2 


/3)(l-« 3/2 ) 


then 


E = -xM 2/3 /2 




where x solves 






= ^x{x) = arccos( 


i — a 2 x) — -\Jl — (1 — a 2 a:) 2 




— arccos(l — aix) + \Jl — (1 — aix) 2 — (t 2 




and 

tB ~ t% x 3/2 






arccos(l — a-ix) — y^l — (1 — a^x) 2 





Elliptic < and at maximum expansion at t? : 
If {t 2 — t 2 ) is close to 

t 2 -h = (a 2 /2) 3/2 [n - arccos(l - 2a) + 2Va - a 2 } 

then a series expansion gives 

E = -xM 2/3 /2 

where x solves 

= ^m{x) ~ 

and 

ts = ti 
w t 2 



93/2 

-2 3 / 2 (2-a 2 x) 1 / 2 + ^-(2-a 2 x) 3 / 2 



+7T — arccos(l — aix) + \Jaix{2 — a±x) — x 3 ^ 2 (t 2 — t±) 



1 



x 



3/2 



arccos 



(1 — a\x) — \Jl — (1 — aix)' 



x~ 3 ' 2 L - 2 3 ' 2 {2 - a 2 xf' 2 + ^(2 - a 2 x) 3 ' 2 



13 



Elliptic E < and recollapsing at t? : 
If 

ti — t\ > (a 2 /2) ' [tt — arccos( a — a 

then 

E = -xM 2/3 /2 

where x solves 



= tpc(x) = 7r — arccos(— 1 + a 2 x) + yl — (1 — a 2 x) 2 

— arccos(l — a\x) + Jl — (1 — aix) 2 — (t 2 — ti)a; 3 ^ 2 



and 



^,3/2 
1 



arccos(l — aix) — J 1 — (1 — aix 



x 



3/2 



7r + arccos(— 1 + a 2 x) + Jl — (1 — a 2 x 



It is easy to adapt the above for the case p 2 > Pi- Clearly any parabolic or hyperbolic 
regions would be collapsing^. 

We conclude this section by stating the result as a theorem: 

Theorem Given any two times t± and t 2 > t±, and any two spherically symmetric density 
profiles < p 2 (M) < pi(M) defined over the same range of M, a L-T model can be found 
that evolves from p x to p 2 in time t 2 — 1\. The inequalities ( ^.ll| )/( p.21[ ) and ( |3.22| )/( p.27| ) 
will tell which class of L-T evolution applies at each M value. The possibilities of shell 
crossings or excessively negative energies are not excluded, and must be separately checked 
for. 



4 Conditions for comoving extrema of density. 

Since we expect the central condensation to propagate outward into the Friedmann back- 
ground, we have to set up the initial conditions so that the edge of the condensation is 
not comoving. For this purpose, it is useful to know the general conditions for comoving 
extrema of density. We shall now consider maxima and minima of p in those domains 
where p is differentiable, and use M as the radial coordinate. 

From ( |2.9D we see that extrema of p will occur at those values of M where 

3 For an expanding and a collapsing hyperbolic region to be contained in the same smooth model, there 
would have to be an elliptic region between them. This is because the p2 < pi region and the p2 > pi 
region must have a point between where pi — pi ■ The only way this can be arranged without causing shell 
crossings is for the elliptic region to be a Kruskal-like neck — see [[Io| 
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(R 3 ),MM — 0. 



(4.1) 



(This is a necessary condition only. Some of the solutions of ( |4.1| ) will be inflection points 
of p(ti,M) rather than extrema, but the whole reasoning below will apply to them,too.) 
For the case E < we find from the second of 



'1 — cos 77)77. 



M 



-2E) 



3/2 



M 



■M 



(t ~ t B ) 



-2Ef' 2 



t B . 



M ■ 



Using this in the first of ( |2.4 ) we find 



(4.2) 



(R 3 ) 



MM 



' M 3 

8E 3 )' 



MM 



1 — cos tjY 



t M 3 " 
M 3 ( 



, M (1 - cosr]) sin 77 



-2E) 



3/2 



sin 2 7] 



8E 3 V 1 — cost? 



+ COS 77 



M 
-2Ef/ 2 



■M 



(t ~ t B ) 



-2E) 



3/2 



M 



■M 



{t ~ t B ) 



M 
-2E) 3 ' 2 



-t 



B,M 



M 3 . 

-3^(1 - cost/) sin 77 



-2E) 3 ' 2 
~M 



At tB,M 



-2E) 3 / 2 
M 

(-2E) 3 / 2 
M 



M 

MM {t — t B ) 



B,M 



B,MM 



0. 



(4.3) 



This equation defines certain values of M, let us call them M = M ex , at which p may have 
extrema. We will verify when they are comoving, i.e. when M ex are independent of time. 

The Jacobian d(M, if) / d(t, M) is nonzero everywhere except those locations where 



V,t=0 



E = 0. 



(4.4) 



Hence, everywhere else M and r\ can be considered to be independent variables. Using the 
second of (|2.4|), we can eliminate (t — t B ) from Q4.3|), and, since the M ex obeying ( |4.3| ) are 
assumed independent of t, what results is an equation in 77 with coefficients depending on 
M. The coefficients of independent functions of 77 all have to vanish. This implies: 



'-2E) 3/2 /M 



,M 



2E) 3 / 2 /M ,MM — t B> M — t Bj MM — 0, 



(4.5) 



(m 3 /^ : 



,MM 



0. 



(4.6) 



all quantities being calculated at M — M ex . Eqs. ( |4.5| ) would have the same form in 
any coordinate system in which M = M(r), but for (|4.6|), the coordinate M is privileged; 
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in other coordinates this equation would look less readable. Note that eqs. (4~o) - (4.6) 
imply that at M = M ex the functions t B , \(-2Efl 2 /Nf\ and (M 3 / E 3 ) agree with their 
Friedmann values (as determined by E 3 / 2 /M and t B at M ex ) up to the second derivatives. 
This means that the local density is at all times the same as that of the Friedmann model 
that matches on there. 

In brief, we have shown that if M = M ex is an extremum of density, the density 
is different iable at M ex , and the extremum is comoving, then (|4.5| ) - ( |4.6j ) are fulfilled. 
Conversely, if p is differentiate at M ex , has an extremum there, and (f4.5|) - ( [4.6D are 
fulfilled, then the extremum will be comoving. 



By the same method it may be verified that in the E(M) > region, where ( |2.6|) apply, 
the conditions for a comoving extremum are again of the form (|4.5| ) - ( |4.6|) , except that 
now the (-2E) in ( |4.5|) is replaced by (2E). 

For the parabolic case, (|2.5|), the condition (R 3 ),mm— reads 



-9 



(2t BM + Mt BtMM ) (t - t B ) - Mt BM 2 =0, (4.7) 



and so the extremum will be comoving when t B> M = ^b,mm = at M = M ex . 

Note that for an arbitrary L-T perturbation inside an exactly Friedmannian exterior, 
the density at the boundary may be discontinuous, as in the Einstein- Straus (12], [L3| config- 
uration. Furthermore, even if it is C 1 at an initial moment, it may develop a discontinuity. 
However, if the above conditions hold at the boundary, they ensure the density is C 1 
through the boundary at all times. 



5 Numerical Example 

5.1 Scales in the background 

The age of the universe is currently believed to be about 14 Gyr. If recombination tem- 
perature is ~ 2700° K, then z = 2700/2.73 ~ 1000 then recombination happened at about 
t r = 3 x 10 5 years when the density was about 10 10 times the present density. 

However, in a k = dust (Friedmann) model, Hq = 65 km/s/Mpc implies to = 2/3Hq = 
10 Gyr, which would put t r at 10 5 yr. 



In an LT model that is close to parabolic even today, eqs (2A) & ( |2.6|) require sinh rj—rj R: 
rj 3 /6 ~ T) — sin?7, so we need rf « 20, say r) < 0.4. For a Friedmann model 

Mp = M r 3 , 2Ep = ±2E r 2 , t BF = (5.1) 

From (|2.4j) & (|2.6|) again the limit on E is 

(2|£ |) 3 / 2 t V 3 0.4 3 , N 

v 1 ou -«-i-< «0.01 5.2 

M 6 6 K ; 
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In a closed model, the maximum in the spatial sections — the 'equator' of the 3-sphere 
- is at r = r m where 2E = — 1, so 

2\E \r 2 m = l =► r m = , 1 (5.3) 

at which point the areal radius today is 

_ ( 9M \ 1/3 ^ 2/3 . . 

^ m °- U(2|^o|) 3 / 2 J ° ( } 

For hyperbolic models, there is no maximum radius, but R m o gives the curvature scale. 
Another way to restrict Eq is to specify that the horizon scale cto be much less than Rmo, 
say t < -Rmo/8; which gives a restriction similar to (|5.2|) . 

A third way to limit E is to specify f2 > 0.03. In a hyperbolic model 

!i= !^ = _^_ (5.5) 

3.R 2 l + cosh(r/) 



Solving this for r\ and using (|2.6|) again gives 

(2|£ |) 3 /% 



< 0.016 (5.6) 



M 

However, within a condensation, the evolution may be nowhere near parabolic. 



5.2 Past null cones, horizons and scales on the CMB sky 



The physical radius of the past null cone in a k = Friedmann dust model with scale 

factor S oc t 2 ^ 3 is t 

L{t) = S [° -dt = 3c (tl /3 t 2 / 3 - t) (5.7) 
Jt S 

so an observed angular scale of 9 on the CMB sky has a physical size at recombination of 

L r = L(t)0 = 3c(tl /3 t 2 r /3 -t r )9 (5.8) 

The present day size of the observed structure — assuming it doesn't collapse — is merely 
scaled up by the ratio of scale factors 

L = L r ^ (5.9) 

The scales would be fairly similar in reasonable k ^ models. 

To determine the condensed structures that correspond to a given present day scale in 
the background, the mass M c of the condensation is divided by the present day density 
Pbfl of the Friedmann background {k = dust) and cube-rooted: 

1/3 



3M C 
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(5.10) 



Conversely, the mass associated with a given scale at a given time is 



M c = v ; - v ; (5.11) 



The particle (causal) horizon at any given time t is 



C = S f* ^dt = 3ct (5.12) 

JO b 



and the visual horizon is 

,t i 



V = S ( -dt = 3c(t-t 1 /h 2/3 ) (5.13) 

Jt r b 

The particle horizon takes no account of inflation, and retains a dust equation of state 
before recombination, so is only included as a rough scale of interest. 



5.3 Scales in the perturbation 

We imagine that present day structures accreted their mass from a background that was 
close to Friedmannian, and therefore the scale of the matter that is destined to end up in 
a present day condensation is fixed by its present day mass. 

The COBE data shows ST/T ~ 10~ 5 on scales of 10°, and the density perturbations 
are Sp/p = 35T/T < 3 x 1CT 5 [|T^]. The power spectrum V{k) = |5k| 2 = where Sp/p = 
Sk<5ke* k ' x is commonly approximated by V(k) = Ake~ ks , where the cut-off scale s is small 
compared to the Hubble scale. This is just V = Ak at longer wavelengths (smaller k) fll6| . 
COBE's measurements had a resolution of ~ 10°, while BOOMERANG'S and MAXIMA'S 
were ~ 0.2°. These angular scales correspond to length scales of 2 Mpc and 50 kpc at the 
time of decoupling, and thus to 2 Gpc and 50 Mpc today. Thus we are only just beginning 
to detect void scale perturbations in the CMB. Although the magnitude of galaxy scale or 
even supercluster scale perturbations, are not yet directly constrained by observations, we 
will retain the figure of ~ 10 -5 . 



The scales associated with present day structures are summarised in the following table. 
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Radius 




Density 


Angle 




today 


Mass 


of sphere 


on CMB 




(kpc) 


(M ) 


(p b ) 


sky (°) 


star 


2 x KT 11 


1 


2 x 10 29 


8 x 10" 7 


globular cluster 


0.1 


10 5 


2 x 10 5 


4 x 10~ 5 


galaxy 


15 


10 11 


6 x 10 4 


4 x lO- 3 


Virgo cluster 


2 000 


2 x 10 13 


5 


0.02 


Virgo supercluster 


15 000 


5 x 10 14 


0.3 


0.06 


Abell cluster (example) 


800 


10 15 


4 000 


0.08 


void 


6 x 10 4 


? 




0.4 


recomb horizon 


280 






1.8 


present horizon 


9.2 x 10 6 






59 


visual horizon 


8.9 x 10 6 






57 


COBE resolution 


1.6 x 10 6 


1.9 x 10 21 


(1) 


10 


BOOM/MAX resolution 


3.1 x 10 4 


1.5 x 10 16 


(1) 


0.2 



Table 1. Approximate scales associated with present day structures. 

Note that the horizons are those in k = Friedmann models without inflation, 
as given in section The masses associated with the resolution scales of 

COBE, MAXIMA & BOOMERANG are obtained by assuming a density equal 
to the parabolic background value pb, as indicated by '(1)' in the density coulmn. 
Useful collections of data can be found at jittp : //www, obspm. f r/messier"7 



tittp : //adc . gsfc.nasa.gov/adc/ sciencedata.html 



and 



tittp : //www . geocities . com/ atlasof theuniverse/supercls .html . 



5.4 Choice of units &; scales 



For the background Friedmann model, we choose the simplest CcLSG, clS its only purpose it 
to get the cosmic timescales & densities approximately right: 



A = 



k = 



p = 



(5.14) 



We will use geometric units such that c — 1 — G, and the remaining scale freedom of 
GR is fixed by choosing units in which the present day mass of the condensation being 
considered is 1. The corresponding geometric length and time units are then: 



M, 



G 



1 



L G = M G G/c 2 



T G = M G G/c 3 



(5.15) 



5.5 The Model 

The principal limitation of the L-T model in the post-recombination era is the absence of 
rotation. However, once rotation has become a significant factor in the collapse process, 
there is already a well defined structure. Later on pressure and viscosity will become 
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important. Our interest is in generating highly condensed structures in a short enough 
timescale, and these factors only come into play once collapse is well underway. Because 
of the lack of rotation etc, all of which tend to delay or halt collapse, we expect our model 
to be rapidly collapsing rather than stationary at the present day. 



We choose to model an Abell cluster: 



MAbell Cluster = 10 15 M Q 
RAbell Cluster = 800 kpC 



From ( |5.15|) and the above table the associated geometric units are 



1 M G 
I L G 

MAbell Cluster 
RAbell Cluster 
t 2 



MAbell Cluster 

48 pc 
156 years 
IM G 
16800 L G 
6.4 x 10 7 T G 



(5.16) 
(5.17) 



(5.11 



5.5.1 Final profile 

At £2 = 10 Gyrs = 6.4 x 10 7 T G , we specify the density profile to be 

p 2 (M)=p M (7000 

which is shown in fig. 4. 



(5.19) 



Fig. 4. goes here 



Now the Friedmann density at t 2 is: 

p b>2 = 1.3 x HT 17 M G /L G = 8 x 10~ 27 kg/m 3 (5.20) 
so the radius in the Friedmann 'background' that contains this mass isQ: 

R F2 = ( 3Ma ^ 1 ciusteX /3 = 26 000 L G (5.21) 

V ^Pb,2 J 

Thus we find 

rM 3 



(R 2 (M)) 3 = / dM> = erfi(4M) (5.22) 

jo 4irp2(M') 224 000V7T p bi2 

as shown in fig. 5, 



For backgrounds with k ^ 0, the radius that contains this mass would be adjusted slightly. 
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Fig. 5 a, b. go here 



and the resulting P2{R) is shown in fig. 6. 



Fig. 6. goes here 



5.5.2 Initial profile 

At t 1 = 100 k years = 10 -5 t 2 = 641 T G we specify the density perturbation to have mass 

M 1 = 10- 2 M Abell cluster (5.23) 

and density enhancement 

3 x 10" 5 p M (5.24) 

for which the chosen profile is: 

/ 1.00003(1 + 100 M)\ , rrtr , 

KM=to{ 1 + 1 ; .003M J (5 ' 25 > 

as plotted in fig. 7. 



Fig. 7. goes here 



The Friedmann density at ti is: 

Pb A = 1.3 x 10~ 7 M G /L% = 8 x 10~ 17 kg/m 3 (5.26) 
and the radius in the 'background' that contains the total mass is: 

R F1 = ( 3MAbM Cluster Y^ = 5700QLg (5>27) 

V ^pb,i J 

The resulting R X {M), 

{R,{M)f = ( M %— dM' = ( M - ln(l + 100 M)) (5.28) 

v v ;; Jo 47rpi(M') 4np bA V 100.003 v ') K ' 

is shown in fig. 8, and pi(R) in fig. 9. 



Figs. 8 & 9. go here 
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5.6 Model Results 



A Maple program was written to generate the formulas and then solve for E(M) and £#(M) 
numerically, as explained in sec. 3. The results are shown in figs. 10. & 11. 



Figs. 10 & 11. go here 



We see that E is of order 10~ 5 which gives a recollapse timescale of 10 7 Tq = 1.7 x 10 9 yr, 
so that the curvature in the condensation is of order Mt 2 /(2E) 3 ^ 2 ~ 0.17. The bang time 
perturbation is of order 2 T G = 300 years, and is quite negligible. 

Stricly speaking, an increasing t B , t B ^ M > 0, creates a shell crossing, but for such a 
slight variation in t B , the shell crossing occurs very early on, long before t\ when the model 
becomes valid. 

The 'velocity' R t — rate of change of the areal radius R — would, in a homogeneous 
model, increase as M 1 / 3 , so plotting i? t /M 1//3 , as in fig. 12, indicates the velocity pertur- 
bation, as a deviation from a constant value. 



Fig. 12. goes here 



In this case, the perturbation is within 3.10 -5 for < M < 0.6, where p 2 is large, but 
increases to 8.10" 4 in the near vacuum region 0.6 < M < 1. This slight excess is most 
likely due to choosing a P2^M) that falls off too fast outside the condensation, requiring 
a too strongly hyperbolic evolution that expands too rapidly. It is likely that this same 
effect requires the slightly increasing t B to keep pi(M) almost flat in these outer regions. 

As a cross-check, these derived functions were used in a separate MATLAB program 
that plots the evolution of a L-T model, given its arbitrary functions. (The appropriate 
form of the evolution equations is given in appendix |B|.) The initial and final density 
profiles were recovered to high accuracy. The resulting density evolution is shown in fig. 
13. 



Fig. 13. goes here 



6 Conclusions 

• We proved that an L-T model can be found to evolve any initial density profile 
on a constant time slice, to any final density profile a given time later. Although 
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it can't be guaranteed the resulting model is free of physical singularities such as 
shell crossings, and the occurrence of a large negative energy (E = —1/2) must 
be handled correctly by changing to M & R decreasing, our numerical experiments 
indicate that realistic choices of the two density profiles and the time difference are 
likely to generate reasonable models. 

• Our numerical example created an Abell cluster in a realistic timescale. It started 
from recombination, with a density perturbation involving a small amount of mass, 
having 5p/p ~ 3.10 -5 . It then 'accreted' most of its final mass. In fact this 'accretion' 
consists of lower expansion rates near the centre, and more rapid expansion rates 
further out. Only at late stages does actual collapse begin at the centre. The initial 
velocity perturbation cannot be chosen if the initial and final densities are chosen. 
It turned out to be Sv/v ~ 3.10 -5 within the future condensation and ~ 8.10 -4 in 
the future vacuum region. The relativiely large value in the outer regions is probably 
due to choosing the final density profile, p 2 to fall off more rapidly than ideal. 

• The preceding two points — the theorem plus the numerical example — demonstrate 
that the L-T model provides a very reasonable description of post-recombination 
structure formation. 

• These two points also indicate that post recombination structure formation in a dust 
universe has an important kinematical component — the initial distribution of ve- 
locities has as much bearing on whether or not a condensation forms and gravity 
magnifies density fluctuations, as the initial density distribution. These initial dis- 
tributions of density and velocity are generated by the functions E(M) and t#(M), 
i.e. coded in the initial conditions. 

• Further numerical examples for structure formation on a variety of scales within the 
L-T model are under investigation. 

• We also obtained the conditions for a density maximum or minimum or shoulder to 
be comoving. Since these are resitrictions on the L-T arbitrary functions, it is evident 
that extrema are in general moving through the fluid, as argued in jT7|], and are not 
comoving. In other words present day density maxima are not likely to be on the 
same worldlines as initial density maxima. 
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A Consistency of the inequalities (|3.21| ) and ( |3.22| ). 



The inequalities G3.21| ) and (|3.22 ) will be consistent if the right-hand side of ( |3.21| ) is smaller 
than the r.h.s. of (|3.22|) in the whole range of a\ and i.e. when 



^ (a 2 3 / 2 - a^ 2 ) < (a 2 /2) 3 / 2 



7r — arccos(l — 2a\/a2) + 2Ja\ja2 — (ai/a2) 2 



(A.l) 



Defining y = ai/a 2 , and recalling that a 2 > a± > by assumption, the above is equivalent 
to 



f(y) > for all < y < 1, where 

4 



Now observe that 



7T - arccos(l - 2y) + 2^y - y 2 - ^(1 - y 2 ) 



/(0)=7T-->0, /(1) = 0, 



(A.2) 



d/ = 2 . 



dy VI - I/ 2 



1-2/-1) <0. 



(A.3) 



Hence, /(y) is monotonically decreasing from /(0) > to and so is positive for all 
< y < 1, which proves ( |A.1| ). 



B Numerical considerations 
B.l Limiting values at M = 0. 

Several of the quantities considered in this paper have the value at the centre of symmetry, 
where M — 0. The variables used in the proof of the theorem, cij and x, have finite limits 
as M — > 0. For numerical programs, these limiting values have to be provided explicitly. 

The values of a« at M = follow very easily. Since a% = R{U,M)/M l / d and R 3 = 
/ M Kp ^ t x } dx, we have, applying the de l'Hospital rule in the third step 



di(0) = lim 



R{t h M) 



1/3 / 



lim — R 3 

M-+0 dM 



lim — 

Af-»o /cp 



lim f ^7 

1/3 / 



3\lV3 



6 



1/3 
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The variable x comes out nonzero automatically when nonzero values of aj(0) are used 
in the program; this follows from the proof of the theorem in sec. |3|. 



B.2 Practical variables 

In practice, it was convenient to define a = a\/a2, Z{ = tj/a2, and y = g^x, and then 
solve the various ip = equations in terms of the variable y, since 02 was in general quite 
large, whereas < ax < 2 in elliptic regions, and ax was at most 200 in a quite strongly 
hyperbolic region. 



B.3 Reconstructing the model evolution 

For reconstructing the evolution of the model, it is convenient to re-write the L-T solutions 

in terms of a; & a: 

Elliptic: 



R = M 1/3 



1/3 (! -cosr/) 



X 



(77 — sin 77) ■ 1 /, 

t-t B = —~ , R = M 1 ' 3 > 



x 



3/2 



\ \ 1 — cos 77 



- 1 



Parabolic or close to it: 



gxi/3 



2800' 



504000 ' 



R = M 1/3 



4[6(« - t B )]-"t ( 1 + |[6(t - t B )]W + ^[6(t - f fl )]*/ 3 - —[(id - /j;)'- 



280' 



3600' 



;b.2) 



Hyperbolic: 



fi = M i/ 3 (cosh77-l) ^ _ ^ = (sinh 77- 77) ^ / _ ) = ^ ] :!( 



x 



3/2 



\ \ cosh 77 — 1 



+ 1 



In all cases the density is 



Ana 2 (a/3 + Ma M ) 
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Fig. 1. The function iPh( x ) for the case E > with ( 3.11 ) fulfilled (middle curve) and with ( 3.11 ) 
not fulfilled (lower curve). The upper curve is the approximating function i/)a(%) from ( |3.12 ). 
The parameters (01,02^2 — *i) are (1,2.5,0.5) for the two upper curves and (1,2.5, 1.4) for the 
lower one. 

1 



0.5 - 



-0.5 




-1.5 

0.2 0.4 0.6 0.8 1 

Fig. 2. The function ipx(%) for the case E < with fljT2lD and (|3T22|) fulfilled (middle curve , 
with (FPU) fulfilled and flOg) not fulfilled (lower curve), and with (|3.21| ) not fulfilled (upper 
curve). The parameters (ai, 02,^2 — h) are (1,2,2) for the middle curve, (1,2,3) for the lower 
curve, and (1,2,0.05) for the upper curve. 
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Fig. 3. The function ipc(%) from Q3.24 ) obeying ( 3.27|) (solid curve) and not obeying it (dashed 
curve). The parameters (ax, a>2, — H) are (1, 2, 9) for the first curve and (1, 2, 0.2) for the other 
one. 
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Fig. 4 The chosen density profile P2(^)/pb2 as a multiple of the 'background' density, for an 
Abell cluster at time ti = 10 Gyr. The axes are in geometric units such that duster = 1) 

as given in (|5.15 ) & ( |5.1£] ). 
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Fig. 5 The areal radius R.2(M) at time t2, that results from the chosen p2(M). The inset shows 
an enlargement of the curve near small M. The axes are in geometric units. 
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Fig. 6 The density profile P2I Pb,2 against areal radius R2. The axes are in geometric units. 
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Fig. 7 The density profile p\{M)/pi,\ chosen for the initial perturbation at time t%. The axes 
are in the geometric units of ( 5.15| ) & (|5.18|) . 




Fig. 8 Areal radius R\{M) at time t\, that is obtained from p\{M). The axes are in geometric 
units. 
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Fig. 9 Density profile pi/pb,i against areal radius R±. The axes are in geometric units. 
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Fig. 10 The L-T energy function E(M) obtained from solving for the L-T model that evolves 
between pi(M) and ^(M). The axes are in geometric units. The symbols "EC", "EX" & 
"H" indicate regions that are respectively elliptic and recollapsing at t2 (EC), elliptic and still 
expanding at ti (EX), and hyperbolic (H). 
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Fig. 11 The L-T bang time function t_e(M) obtained from solving for the L-T model that 
evolves between pi(M) and p2(M). The axes are in geometric units. 
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Fig. 12 The velocity perturbation R/M 1 / 3 at time t\. A constant value would indicate no 
perturbation. The axes are in geometric units. 
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Fig. 13 The evolution of p(t, M) for the derived L-T model. The axes are in the geometric units 
of ( gig) & (|5l8|) . In the range < M < 0.795 the evolution is elliptic and already recollapsing 
at time £2, in 0.795 < M < 0.865 it is elliptic but still expanding at £2, and for M > 0.865 it 
is hyperbolic. In practice, recollapse would be halted at some point by the effects of pressure, 
rotation, etc. The initial and final density profiles calculated at times ti and t<i coincide with 
those originally chosen and shown in figures 7 and 4 respectively. 
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